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1  RESEARCH  PROGRESS  AND  ACTIVITIES 
REPORT 

During  the  last  eighteen  months  of  the  grant  our  work  concentrated  on  two 
main  research  topics: 

(i)  “Parallel  and  Sequential  Iterative  Methods  for  Linear  and  Non¬ 

linear  Systems”.  Much  of  the  work  on  this  topic  concentrated  on  the 
convergence  and  rate  of  convergence  of  parallel  asynchronized  methods 
for  .'ol'dng  linear  systems  mising,  on  the  one  hand,  from  the  numer¬ 
ical  solution  to  partial  differential  equations  and,  on  the  other  hand, 
from  least  squares  solution  to  rectangular  systems  which  arise  in  ap¬ 
plication  such  as  image  reconstruction  from  incomplete  tomographical 
data.  The  mathematics  behind  the  analysis  of  these  two  applications 
of  the  asynchronized  parallel  methods  is  quite  different. 

Recently  we  have  been  able  to  extend  our  convergence  results  to  asyn¬ 
chronized  methods  for  solving  nonlinear  systems.  One  application  now 
consists  of  tomographic  reconstruction  from  incomplete  data  where  the 
image  is  constrained  to  lie  in  a  bounded  convex  set  such  as  an  n  di¬ 
mensional  box. 

(ii)  “Reachability  Problems  for  Dynamical  Systems”.  Here  we  con¬ 

centrated  on  developing  numerical  methods  to  test  whether  the  trajec¬ 
tory  of  a  linear  differential  system  emanating  from  a  given  initial  state 
becomes  from  some  point  onwards  nonnegative.  We  particularly  char¬ 
acterized  when  such  intial  states  are  symbioisis  points,  meaning  that 
from  some  peunt  onwards  all  populations  become  nondecreasing. 

As  a  by-product  of  the  work  on  these  topics  we  also  had  to  sedve  various 
theoretical  problems  which  can  be  described  under  the  heading: 


(iii)  “Problems  in  Nonnegative  Matrices  and  their  Applications”. 

We  shall  describe  the  main  results  which  were  obtained  on  these  topics  in 
the  next  section  of  this  report.  We  strongly  believe  that  an  examination  of 
the  restdts  which  were  achieved  over  the  life  of  the  grant  shows  that  many  of 
the  goals  which  were  suggested  in  the  intial  1987  proposal  and  in  the  annual 
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and  progress  reports  which  have  been  submitted  since,  have  been  realized. 
Quite  a  few  of  the  questions  that  have  been  raised  have  been  answered,  but 
not  always  with  the  solution  that  was  conjectured. 

Since  the  beginning  of  the  work  on  the  proposal,  15  papers  which  sum¬ 
marize  our  results  on  the  above  three  topics  have  been  submitted  for  publi¬ 
cation.  Their  titles  are  as  follows; 

On  parallel  and  sequential  iterative  methods 

1)  (with  P.  J.  Kavanagh)  “Consistency  and  convergence  of  the  parallel  mul¬ 

tisplitting  method  for  singular  M-matrices,”  SIAM  J.  Matrix  Analysis 
Appl.,  10(1989),  pp.  210-218. 

2)  (with  L.  Eisner  and  I.  Koltracht)  “On  the  convergence  of  asynchronous 

paracontractions  with  application  to  tomographic  reconstruction  from 
incomplete  data,”  Lin.  Alg.  Appl.,  130(1990),  pp.65-82, 

3)  (with  A.  Hadjidimos)  “Convergence  domains  of  the  SSOR  method  for 

generalized  consistently  ordered  matrices”,  J.  Comp.  Appl.  Math., 
33(1990),  pp.  35-52, 

4)  (with  M.  Hanke)  “Preconditioning  and  splittings  for  rectangular  sys¬ 

tems”,  Numer.  Math.,  57(1990),  pp.8^95. 

5)  (with  E.  Eisner  and  B.  Vemmer)  “The  effect  of  the  number  of  processors 

on  the  convergence  of  the  parallel  block  Jacobi  method”,  Lin.  Alg. 
Appl.,  154-156(1991),  pp.3U-330. 

6)  (with  M.  Hanke  and  W.  Niethammer)  “On  the  SOR  method  for  sym¬ 

metric  positive  definite  systems”,  Lin.  Alg.  Appl.,  154-156(1991), 
pp.457-472. 

7)  (with  L.  Eisner)  ”Monotonic  sequences  and  rates  of  convergence  of  asyn¬ 

chronized  iterative  methods”,  submitted  to  Numer.  Math. 

8)  (with  L.  Eisner  and  I.  Koltracht)  “Convergence  of  sequential  and  asyn¬ 

chronous  nonlinear  paracontractions”,  submitted  to  Numer.  Math. 

On  the  reachability  problem 

9)  (with  R.  J.  Stern)  “  Discrete  approximations  to  reachability  cones  of 

linear  differential  equations,”  Lin.  Alg.  Appl.,  120(1989),  pp.  65-79. 
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10)  (%ith  R.  J.  Stern  and  M.  Tsatsomeros)  “The  reachability  cones  of  essen¬ 

tially  nonnegative  matrices”,  Lin.  Multilin.  Alg.,  28(1991),  pp.213- 
224. 

11)  (with  M.  Tsatsomeros)  “Symbiosis  points  for  linear  differntial  systems”, 

Lin.  Multilin.  Alg.,  30(1991),  pp.49-59. 

On  nonnegative  matrices  and  applications 

12)  (with  R.  E.  Hartwig  and  N.  J.  Rose)  “  An  algebraic-analytic  approach 

to  nonnegative  basis,  ”  Lin.  Alg.  Appl.,  133(1990),  pp.77-88. 

13)  (with  I.  Koltracht)  “On  the  inverse  M-matrix  problem  for  real  symmet¬ 

ric  positive  definite  Toeplitz  matrices”,  SIAM  J.  Matrix  Anal.  Appl., 
12(1991),  pp.310-320. 

14)  (with  H.  J.  Werner)  “Nonnegative  group  inverses”,  Lin.  Alg.  Appl., 

151(1991),  pp.85-96. 

15)  (with  H.  Schneider)  “Principal  components  of  minus  M-matrices”,  sub¬ 

mitted  to  Lin.  Multilin.  Alg. 


During  the  first  half  of  the  period  in  which  this  grant  has  been  in  effect 
we  have  also  co-authored  a  book  in  connection  with  the  second  research 
topic  listed  above.  The  title  of  the  book  is  “  Nonnegative  Matrices  in  Dy¬ 
namic  Systems”.  Its  other  authors  are  A.  Berman  and  R.  J.  Stern  and  it 
was  published  in  the  Series  in  Pure  and  Applied  Mathematics,  Wiley  Inter- 
science,  New  York,  1989. 

Since  the  start  of  the  grant  we  have  presented  results  from  our  work  and 
acknowledged  the  support  of  the  Air  Force  in  conferences  and  coUoquia  as 
follow:  “Special  Session  on  Modern  Trends  in  Matrix  Theory  and  its  Ap¬ 
plications”,  AMS  Annual  Meeting,  January,  1988;  “Oberwolfach  Meeting 
on  Numerical  Algebra  and  Parallel  Computations”,  Oberwolfach,  February 
1988;  “Workshop  on  Iterative  Solutions  to  Singular  Systems”,  Univ.  of  Karl¬ 
sruhe,  Karlsruhe,  West  Germany,  March  1988;  ”SIAM  3rd  Conference  on 
Applied  Linear  Algebra”,  Madison,  May  1988;  “International  Symposium 
on  Computational  Applied  Math.”,  Leuven,  Belgium,  July  1988;  “NATO 
Advanced  Study  Institute  on  Numerical  Linear  Algebra,  Digital  Signal  Pro¬ 
cessing  and  Parallel  Algorithms”,  Leuven,  Bel^um,  August  1988;  “Confer¬ 
ence  on  Iterative  Methods  for  Large  Linear  Systems  (dedicated  to  David 
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M.  Young)”,  University  of  Texas,  Austin,  October  1988;  “A  Conference  on 
Approximation  Theory  and  Numerical  Linear  Algebra  (dedicated  to  R.  S. 
Varga)”,  Kent  State  University,  Kent,  Ohio,  March  1989;  Northern  Illinois 
University  Conference  on  “Linear  Algebra,  Numerical  Linear  Algebra  and 
Applications”,  DeKalb,  Illinois,  April  1989;  “The  Householder  Symposium 
XI  on  Numerical  Linear  Algebra”,  Tylosand,  Sweden,  June  1990;  Meeting 
on  “Numerical  Linear  Algebra”,  Oberwolfach,  Germany,  April  1991;  “Haifa 
VII  Matrix  Theory  Conference”,  Haifa,  Israel,  June  1991. 

The  referencing  within  the  report  is  as  follows.  Papers  cited,  but  not 
co-authored  by  the  P.I.,  are  referenced  by  numbers  in  the  text  and  the  key 
is  given  after  section  3.  References  to  papers  co-authored  by  the  P.I.  are 
numbered  by  [Nxx],  where  xx  refers  to  the  paper  number  in  the  P.I.’s  vitae 
which  is  attached  at  the  end  of  this  report. 
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2  ASYNCHRONOUS  ALGORITHMS  FOR  LARGE 
LINEAR  SYSTEMS 


In  this  section  we  describe  our  research  concerning  the  convergence  and  the 
acceleration  of  convergence  of  a  certain  model  of  a  parallel  chaotic  (also 
known  as  asynchronized)  iteration  scheme.  Some  of  the  problems  we  en¬ 
countered  come  from  the  fact  that  we  tried  to  apply  the  same  asynchronized 
model  to  linear  systems  whose  coefficient  matrices  arises  in  different  appli¬ 
cations  and  also  to  nonlinear  systems.  This  means  that  for  each  application 
we  had  to  find  the  inherent  mathematical  properties  which  make  the  conver¬ 
gence  of  the  algorithm  and  its  acceleration  possible.  Each  type  of  system, 
in  turn,  gives  rise  to  different  problems  in  the  actual  implementation  of  the 
algorithm. 

The  chaotic  iteration  method  which  we  have  in  mind  has  the  fcilowing 
form:  We  are  given 

i)  m  linear  or  nonlinear  operators  Si, ... ,  Sm- 

ii)  A  computation  cycle,  namely,  a  fixed  time  period  T  >  0,  and  a  reg¬ 

ulated  sequence  integers  on  m,  that  is  a  sequence  of  integers  {»i}^i  with 
1  £  and  such  that 

{l,2,...,m}C  .,»,+r-i},Vj>l.  (2.1) 

iii)  m  nonnegative  diagonal  matrices  Ei,  I  =  1,. . .  ,m,  whose  sum  is  the 
identity  matrix.  (They  are  sometimes  called  weighting  or  masking  matrices.) 

iv)  A  parallel  machine  with  k  processors  and  a  host  node. 

We  perform  the  iteration: 

=  (7  -  4-  j  =  1,2,....  (2.2) 

Our  model  works  as  follows:  At  time  j  a  processor,  call  it  for  now  the 
subject  processor,  which  has  just  completed  a  previous  task  is  assigned 
the  task  specified  by  ij,  namely,  by  the  operators  and  This  means 
that  it  begins  to  calculate  the  vector  u  =  Ei^Bi^x^^K  Note  that  only  the 
entries  of  u  corresponding  to  the  nonzero  diagonal  entries  in  Ei^  need  be 
computed.  The  number  r^  -  1,  <  T,  then  represents  the  number  of 
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similar  tasks  completed  by  other  processors  before  the  subject  processor 
completes  its  present  computation.  When  this  computation  is  done,  the 
sum  u  +  (/—  Ei^  is  formed  by  the  host  processor  and  the 

subject  processor  is  then  assigned  task  ij+tj-  The  existence  of  a  computa¬ 
tion  cycle  r  >  0  as  given  by  (2.1)  means  that  there  exists  a  time  period  T 
such  that  in  every  T  successive  iteration,  the  global  approximation  to  the 
solution  is  corrected  at  least  once  by  each  of  the  operators  Bi,... ,  Bm- 

We  have  proved  the  convergence  of  (2.2)  when  the  chaotic  process  is  used 
to  solve  iteratively  two  quite  different  types  of  linear  systems 

Ax  =  b.  (2.3) 

The  first  type  usually  arises  from  a  finite  differences  approximation  to  second 
order  partial  differential  equations  subject  to  boundary  value  conditions. 
There  A  is  frequently  a  monotone  matrix,  meaning  it  is  nonsingular  with 
A~^  >  0.  The  operators  5/,  I  =  are  then  iteration  matrices 

induced  by  m  weak  regular  splittings  of  the  matrix  A,  that  is,  by  m  splittings 

A  =  Ml  -  Ni,  (2.4) 

satisfying  Mi  is  invertible  with 

M-‘  >0  and  Bi  =  Mf^Ni  >  0.  (2.5) 

The  second  type  of  linear  systems  to  which  we  have  applied  (2.2)  are  rect¬ 
angular  systems  (2.3)  which  arise  in  image  reconstruction  from  incomplete 
data  as,  for  example,  in  well-to-well  tomography  used  in  geophysics.  Previ¬ 
ously  the  cyclic  Algebraic  Reconstruction  Technique  (ART),  which  itself  is  a 
generalization  of  the  Kaczmarz  projection  method,  has  been  applied  to  find 
the  least  squares  solution  of  minimal  norm  to  such  systems.  The  cyclic  ART 
method,  which  is  also  closely  related  to  the  successive  overrelaxation  (SOR) 
method  (see  Koltracht  and  Lancaster  (1]  and  Hanke  and  Niethaminer  [zj), 
is  a  sequential  method  where  we  apply  in  a  cyclic  order  the  m  operators 
Bi,  each  of  which  corresponds  to  an  orthogonal  projection  onto  a  subspace 
spanned  by  a  row  or  a  group  of  rows  of  the  coefficient  matrix  A  with  each 
row  of  the  A  appearing  in  at  least  one  of  the  groups  (thus  overlapping  is 
allowed).  The  proof  that  the  cyclic  ART-SOR  can  be  parallelized  according 
to  (2.2)  is  more  intricate  than  in  the  case  of  chaotic  iteration  for  s(Jving 
monotone  system.  It  requires  certain  norm  considerations  and  restrictions 
on  the  weighting  diagonal  matrices  which  are  not  necessary  in  the  case  of 
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chaotic  iterations  for  monotone  systems. 

For  both  types  of  linear  systems  mentioned  above  we  have  found  various 
proofs  for  the  convergence  of  (2.2)  all  of  which  involve  the  embedding  of 
the  process  as  a  sequential  iteration  process  which  takes  place  in  higher  di- 
mensionsal  space.  In  one  of  the  types  of  embedding,  we  iterate  sequentially 
in  the  fcn-dimensional  space  (recall  k  is  the  m  mber  of  processors  and  n  is 
the  dimension  of  the  iterates  in  (2.2))  and  produce  a  sequence  of  iteration 
vectors  =  1,2, _  The  idea  now  is  to  prove  that  the  k  subvectors  of 

G  i?"*  each  has  a  limit  point,  as  j  — »  oo,  equal  to  the  stJution  of  (2.3). 
Although  we  have  not  used  this  name  formally  in  any  of  our  reports,  we 
like  to  refer  informally  to  this  approach  to  the  proof  as  the  “logbook” 
approach.  This  is  because  for  any  processor,  say  the  jz-th,  we  think  of  the 
sequence  of  n-vectors  which  we  can  form  from  the  t/-th  subvectors  of  the 
z^^^  as  keeping  a  “logbook”  of  the  local  approximations  in  the  i/-th  proces¬ 
sor  at  each  time  step  of  the  global  iteration.  Thus  much  of  the  time  only 
the  index  of  the  iteration  in  the  processor  is  advanced,  but  the  actual  value 
of  the  local  approximation  is  unchanged.  It  only  changes  when  the  iz-th 
processor  updates  to  and  downdates  from  the  host  processor. 

The  second  type  of  embedding  we  have  used  to  prove  the  convergence  of 
(2.2)  is  by  blowing  up  the  process  in  the  n-dimensional  space  to  a  sequential 
process  in  the  nT-dimensional  space,  where  T  is  the  computation  cycle.  This 
is  achieved  by  looking  at  the  iteration 


where 


C 


rj  blocks 

'/(/-£.,)  0  ...  0 

/  .  0  0 
0  / 

0 

0  .  0  I 


V  0  .  0 


(2.6) 


/  0/ 
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and  where 


y 


(») 


Above,  under  the  appropriate  assumptions  on  A  and  the  fl/’s,  a  contrac¬ 
tion  occurs  at  least  every  time  a  succession  of  2T  iterations  has  been  applied. 
In  otherwords,  in  an  appropriately  induced  norm  on  HCfc+^T-i  •  •  •  C’ltH  < 


1,  Vit  >  1. 


One  of  the  goals  of  our  research  was  to  obtain  a  better  understanding 
of  the  meaning  and  the  interpretation  of  the  two  embeddings.  First,  we 
can  view  rj  as  che  time  elapsed  (=  time-  lag)  between  two  updates  of  the 
global  approximation  by  the  operator  ij.  Thus  if  rj  is  large,  then  the  sub¬ 
ject  processor  which  received  the  global  approximation  is,  by  the  time 
it  has  finished  computing  u,  using  a  correction  based  on  a  relatively  old 
approximation  to  update  the  current  global  approximation  in  the  host 
node.  Using  special  nested  subcones  of  monotonic  vectors  in  ,  the  cone 
of  nonnegative  vectors  in  the  i2”^-climensional  space,  we  were  able  to  prove 
in  [N61)  the  fcAlowing  result  on  the  rate  of  converegence  of  the  model  pven 
in  (2.2). 


Theorem  1  Let  and  be  two  sequences  of  time-laggs  such 

that 

1  <  <rj  <T,  j  =  1,2 .  (2.7) 


If 


+  l  ^  ^j"bl,  j  1(2,..., 


(2.8) 


then 


sup  limsup  IjCj  . .  .Ciy  -  <  sup  limsup  ||C;  . . . Ci y  -  (2.9) 

J— oo  yef^^  J— oo 
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where 


Sj  blocks 


and  where  ^  =  (i^, . . .  G  with  x  being  the  solution  to  (2.3). 

Let  us  refer  to  the  iterative  processes  induced  by  the  sequences  of  time- 
laggs  and  {rj}^j  as  the  more  frequently  updtaing  process  and 

the  more  infrequently  updating  process,  respectively.  What  the  above 
result  says  is  this:  When  5jo+i  >  ^jo  +  ^  some  jo  >  1,  then  the  more  fre¬ 
quently  updating  process  uses  an  older  approximation  to  compute  the  jo-th 
iteration  than  the  approximation  it  has  used  in  computing  the  immediately 
proceeding  iteration.  Therefore  condition  (2.8)  means  that  when  the  more 
frequently  updating  process  never  “suddenly”  uses  an  cJder  approximation 
in  computing  some  iterate  than  the  approximation  it  has  used  in  computing 
the  previous  iterate,  then  the  rate  of  convergence  of  the  more  frequently 
updating  iteration  is  more  favorable  than  the  rate  of  '’onvergence  of  the 
more  infrequently  updating.  We  have  shown  by  means  of  examples  that  if 
condition  (2.7)  holds,  but  condition  (2.8)  does  not,  then  the  result  of  the 
theorem  is  not  true. 

We  have  carried  out  many  numerical  experiments  in  connection  with 
Theorem  1  and  several  such  are  given  in  [N57].  There  we  considered  the 
special  case  when: 

Sj  =  k  =  const,  and  rj  =  k'  =  const.,  Vj  >  1 

and  with 

it  <  it'. 
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Such  a  delay  will  occur  when  the  work  among  the  processors  is  equally 
distributed  in  which  case  the  constant  delay  is  just  the  number  of  processors 
minus  1.  In  the  following  table,  for  a  typical  80  X  80  diagonally  dominant 
matrix,  we  shew  how  the  increasing  the  number  of  processors  (viz.  increasing 
the  delay)  effects  the  number  of  iteration  which  are  necessary  to  reach  a 
0ven  accuracy  of  10“®.  In  the  table  k  is  the  number  of  processors  used  by 
the  machine  and  J  =  J{k)  is  the  number  of  iterations  required  to  achieve 
a  prescribed  accuracy  to  the  solution.  The  last  column  is  the  ratio  of  the 
number  of  iterations  to  the  number  of  processors,  indicating,  roughly,  the 
number  of  iterations  which  each  of  the  processors  would  have  to  execute  in 
parallel  if  communication  overheads  are  reasonable: 


k 

J 

J/k 

k 

J 

J/k 

1 

101 

101.0 

11 

291 

26.5 

2 

134 

67.0 

12 

324 

27.0 

3 

162 

54.0 

13 

352 

27.1 

4 

194 

48.5 

14 

384 

27.4 

5 

204 

40.8 

15 

394 

26.3 

6 

212 

35.3 

16 

402 

25.1 

7 

244 

34.9 

17 

434 

25.5 

8 

261 

32.6 

18 

454 

25.3 

9 

282 

31.3 

19 

474 

24.9 

10 

291 

29.1 

20 

484 

24.2 

To  give  a  better  illustration  of  this  table  of  iterates  let  us  graph  k  versus 
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100  - 

80  - 

60- 

40- 
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The  graph  clearly  points  to  a  conjecture  that  we  have  that  aa  the  number 
of  processors  increase  the  speed-up  tends  to  a  linear  constant  and  no  gain 
is  achieved  by  increasing  heavily  the  number  of  processors. 

For  rectangular  systems  of  equations  we  have  considered  chaotic  itera¬ 
tions  for  computing  least  squares  solutions.  As  mentioned  erlier,  the  math¬ 
ematics  that  is  needed  to  demonstarte  that  such  iterations  converge  is  quite 
different  than  for  monotone  systems.  We-attach  to  this  report  a  recent  paper 
in  which  we  show  that  chaotic  methods  of  the  form  (2.2)  can  be  applied  also 
to  nonlinear  systems  of  equations.  This  allows  us  also  to  consider  applica¬ 
tions  to  finding  linear  least  squares  solutions  lying  in  some  closed  convex  set 
which  represents  a  nonlinear  constraint  on  the  solution.  Such  an  application 
arises  in  computed  tomography  from  incomplete  data. 
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3  The  Determination  of  Nonnegative  Solutions 
to  Linear  Systems  of  Differential  Equations  by 
Finite  Differences  Methods 

Consider  the  system  of  linear  differential  equations 

i  =  Axy  (3-1) 

where  A  =  (a,j)  is  an  n  X  n  real  matrix.  In  many  engneering,  biological, 
and  other  applications  the  vector  x{t)  represents  the  state  of  a  system  at 
time  t  and  its  components  frequently  represent  the  sizes  of  populations  or 
species  at  time  t.  In  some  applications  (see  Luenberger  [1979])  the  matrix 
A  is  essentially  nonnegative,  that  is,  >  0,  *  51^  j.  Such  a  constraint  en¬ 
sures  that  trajectories  which  emanate  from  nonnegative  initial  states  remain 
nonnegative.  Moreover  interest  centers  on  trajectories  which  eventually  be¬ 
come  and  remain  nonnegative  or  trajectories  whose  velocities  (derivatives) 
become  and  remain  nonnegative.  If  the  latter  condition  holds,  then, 
in  time,  the  system  reaches  a  state  from  which  every  species  will 
not  decrease  in  size  thereafter.  We  call  initial  points  whose  trajectories 
and  their  velocities  eventually  become  and  remain  nonnegative  it  symbio¬ 
sis  points. 

From  now  on  we  shall  suppose  that  A  is  essentially  nonnegative.  Denote 
by  Xa{R+)  the  set  of  all  pcsnts  in  IV'  such  that  the  trajectories  emanat¬ 
ing  from  these  points  become  and,  due  to  the  essential  nonnegativity  of  A, 
remain  nonnegative.  In  [N33]  we  showed  that  Xa(R^)  is  a  convex  cone 
which,  however,  need  not  be  closed  or  pointed.  In  a  sequence  of  papers 
[N31],  [N33],  and  [N37]  we  gave  formulas  for  the  closure  of  X/i(ii!].)  under 
various  further  assumptions  on  A  such  as  diagonalizability,  real  spectrum, 
etc.  These  formulas  were  very  difficult  to  apply  for  two  reasons:  (i)  they 
were  too  complicated  as  they  involved  the  intersections  of  the  eigenspaces 
of  A  with  various  projections  of  the  nonnegative  orthant,  and  (ii)  as  only 
the  closure  of  the  reachability  cone  was  determined,  there  were  further  com¬ 
plications  in  applying  the  formulas  to  determine  whether  a  given  boundary 
point  of  Xa(R^)  is  also  a  reachability  point. 

Because  of  the  difficulties  we  described  above  we  thought  of  the  possi¬ 
bility  of  applying  finite  differences  methods  in  order  to  determine  whether  a 
given  pcMnt  xq  G  i?"  is  a  reachability  point.  In  the  simplest  finite  differences 
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schemes,  the  so  called  Cauchy-Euler  method,  we  approximate  the  solution 
at  times  k  =  1,2, .. .  from  the  quotients 

(3.2) 

where  h  is  the  time-step  used  by  the  method.  After  some  rearrangement 
we  obtain  from  (3.2)  the  discrete  trajectory  o/pomt5  emanating  from  xo  =  xo 
given  by 

xk  =  {I^hAfxo,  /:  =  1,2,....  (3.3) 

We  see  that  discretization  schemes  for  systems  of  ordinary  differential  equa¬ 
tions  thus  resemble  error  analysis  for  iterative  solutions  to  linear  systems  of 
equations  in  the  sense  that  both  procedures  involve  powering-up  matrices. 
It  is  therefore  of  no  surprise  that  underlying  both  are  basic  features  and 
problems  of  the  power  method  for  determining  eigenvalues  and  invariant 
subspaces. 

Observe  that  if  the  time  step  h  is  small  enough  to  make  the  matrix 
I  -k-hA  nonnegative,  then  if  xq  €  is  a  point  for  which  there  exists  an  an 
exponent  ko  such  that  Xko  is  nonnegative,  then  all  subsequent  points  in  the 
trajectory  emanating  from  xq  will  remain  nonnegative.  One  result  that  we 
completed  proving  during  the  course  of  this  grant  represents  a  considerable 
improvement  over  results  which  we  obtained  previously  in  [N48].  It  is  the 
following: 

Theorem  2  ([N54])  Let  A  be  an  n  x  n  essentially  nonnegative  matrix  and 
consider  the  linear  differential  system  (3.1).  Let 

/i(A)  :=  sup{h  >  0  I  / -I- /lA  >  0}.  (3.4) 

Then  xq  E  Ayi(iZ"  )  if  and  only  if  for  any  0  <  h  <  h{A)  there  exists  an  index 
ko  such  that  the  discrete  trajectory  of  points  (3.3)  generated  from  xq  =  xq 
satisfies  that 

Xk  >  0?  yk  >  ko. 

Notice  that  h{A)  can  be  very  large,  it  is  -boo  if  A  is  nonnegative,  but 
in  any  case  it  depends  only  in  the  size  of  the  diagonal  entries  of  A  and  is 
not  “infinitesimal”.  What  the  result  means  is  this:  “regardless  of  the  ex¬ 
tent  to  which  the  continuous  and  discrete  trajectories  diverge  from 
each  other,  one  becomes  nonnegative  if  and  only  if  the  other  one 
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does,  provided  only  that  the  time— step  h  satisfies  0  <  h  <  h{AY. 
This  is  a  qualitative  well  as  a  numerical  statement  about  the  behavior  of 
the  solutions  to  systems  of  ordinary  differential  equations  whose  coefficient 
matrix  is  essentially  nonnegative.  The  proof  of  the  above  theorem  is  quite 
involved  and  is  the  sul_  t  of  the  manuscript  “Reachability  cones  of  essen¬ 
tially  nonnegative  matrices”  which  has  just  been  accepted  for  publication  in 
the  journal  of  Linear  and  Multilinear  Algebra.  We  mention  that  quite  an 
important  tool  which  was  used  in  the  proof  of  the  theorem  is  taken  from 
an  earlier  paper  [N52]  in  which  we  considered  an  analytic  approach  to  the 
question  of  existence  of  a  nonnegative  basis  for  the  eigenspace  of  a  nonneg¬ 
ative  matrix  corresponding  to  its  Perron  root. 

The  characterization  of  symbiosis  points  is  done  in  [N59].  Decompose  a 
pcant  u  G  Xa{R^)  into 

V  =  t;+  — 

where  is  the  join  of  all  eigenspaces  of  A  corresponding  to  eigenvalues 
with  a  nonnegative  real  part  and  u_  is  in  the  join  of  all  eigenspaces  of  A 
corresponding  to  eigenvalues  with  a  negative  real  part.  Thus  t;_  is  in  the 
stability  part  of  the  space  since  limt-»oo  =  0-  v  we  define  the 

invariant  set  of  components  of  i;  as  the  set 

I{v)  =  {1  <  t  <  n  :  =  (v+),-,  Vt  >  0  }.  (3.5) 

Thus  I(u)  consists  of  the  indices  of  the  components  of  the  vector  which 
remain  invariant  throughout  the  entire  trajectory  emanating  from  We 
prove  the  following  characterization: 

Theorem  3  ([N59])  Let  Abeannxn  essentially  nonnegative  matrix.  Then 
a  vector  v  £  )  is  a  symbiosis  point  for  (3.1)  if  and  only  if  there  exists 

a  sufficiently  large  time  to  such  that 

j  G  T{v)  =>  0  <  (e*^v_)j  j,  0,  Vt  >  to,  (3-6) 

where  T{v)  is  given  in  (3.5).  Furthermore,  if  v  £  Xa{R^)  is  a  symbiosis 
point,  then  for  any  j  £  X{v),  (u+))  =  0  if  and  only  if  (e‘"^t))j  =  0  for  all 
t  >  0. 

We  have  two  comments.  First,  in  the  spirit  of  Theorem  2  we  can  also 
characterize  a  symbiosis  point  u  £  XA{Lt^)  in  terms  of  the  nondecreasing- 
ness  of  finite  differences  sequences  generated  from  v  (similar  to  the  way 
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in  which  the  sequence  in  (3.3)  is  generated  from  xo)  where  the  step  size 
h  satisfies  (3.4).  Second  in  the  case  when  A  is  weakly  stable,  meaning 
its  eigenvalue  with  the  largest  real  part  is  the  origin,  then  symbiosis  points 
admit  a  matrix-combinatorial  structure  in  the  sense  that  it  is  possible  to 
determine  apriori  which  indices  1  <  i  <  n  lie  in  I(u)  from  a  certain  block 
directed  graph  of  the  matrix  A.  Both  of  these  issues  are  addressed  in  [N59]. 
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